Pb 2

Pb 2#

try:
    import dolfinx
    print("dolfinx importé avec succès")
except ImportError as e:
    print(f"Erreur lors de l'importation de dolfinx : {e}")
    import sys
    print(f"Python path : {sys.path}")

dolfinx importé avec succès

Problème torsion d’arbres#

Librairies#

# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue
import ipywidgets as widgets
from IPython.display import display
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')

Géométrie#

Définition de la géométrie et du maillage#

!{thebe-button}
# Initialisation du plotteur PyVista en dehors de la fonction
p = pv.Plotter(notebook=True)
p.set_background("grey")

# Fonction pour créer et visualiser la géométrie
def create_ellipse(rho_x, rho_y, h=1.0, lc=0.02):
    global domain
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("ellipse")

    # Création de l'ellipse
    ellipse = gmsh.model.occ.addEllipse(0, 0, 0, rho_x, rho_y)
    curve_loop = gmsh.model.occ.addCurveLoop([ellipse])
    surface = gmsh.model.occ.addPlaneSurface([curve_loop])

    # Extrusion pour obtenir un cylindre elliptique
    gmsh.model.occ.extrude([(2, surface)], 0, 0, h)
    gmsh.model.occ.synchronize()

    # Configuration du maillage
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
    gmsh.model.mesh.generate(3)

    # Sauvegarde et conversion du maillage
    gmsh.write("ellipse.msh")
    msh = meshio.read("ellipse.msh")
    meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))
    gmsh.finalize()

    # Lecture du maillage avec FEniCSx
    with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
        domain = xdmf_file.read_mesh(name="Grid")
        domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
    
    # Conversion du maillage pour la visualisation avec PyVista
    u_topology, u_cell_types, u_geometry = plot.vtk_mesh(domain)
    u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)

    # Effacer la scène précédente et ajouter le nouveau maillage
    p.clear()
    # Ajout du maillage à la scène de visualisation
    p.add_mesh(u_grid, 
               show_edges=True,  # Affiche les arêtes du maillage
               scalar_bar_args={
                   "title": "u",  # Titre de la barre de couleur
                   "title_font_size": 24,
                   "label_font_size": 22,
                   "shadow": True,
                   "italic": True,
                   "font_family": "arial",
                   "vertical": False  # Orientation horizontale de la barre de couleur
               })

    # Ajout d'un titre à la visualisation
    p.add_text("Cylindre Mesh", font_size=24, color="black", position="upper_edge")

    # Ajout des limites de la boîte englobante
    p.show_bounds(color="black")

    # Ajout des axes de coordonnées
    p.add_axes(color="black")

    # Définition de la couleur de fond
    p.set_background("grey")
    
    # Mise à jour de la scène
    p.show()
    
    return rho_x, rho_y, h
    
# Widgets interactifs pour a, b, et h
rho_x_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-grand axe')
rho_y_slider = widgets.FloatSlider(value=0.2, min=0.1, max=1.0, step=0.01, description='Demi-petit axe')
h_slider = widgets.FloatSlider(value=1.0, min=0.5, max=3.0, step=0.1, description='Hauteur')

# Interface utilisateur et fonction interactive
ui = widgets.VBox([rho_x_slider, rho_y_slider, h_slider])
out = widgets.interactive_output(create_ellipse, {'rho_x': rho_x_slider, 'rho_y': rho_y_slider, 'h': h_slider})

display(ui, out)

# Fonction pour récupérer les valeurs actuelles des curseurs
def get_slider_values():
    current_rho_x = rho_x_slider.value
    current_rho_y = rho_y_slider.value
    current_h = h_slider.value
    return current_rho_x, current_rho_y, current_h

# Exemple d'utilisation : récupération des valeurs
rho_x, rho_y, h = get_slider_values()
rho = np.array([rho_x, rho_y])  # Demi-axes dans les directions x et y

zsh:1: command not found: thebe-button

Configuration du problème#

Définition de l’espace fonctionnel#

# Définition de l'espace fonctionnel pour le problème
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

Définition des frontière du domaine#

Visualisation des frontière du domaine#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites#

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc0_S0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc0_Sh = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc0_Sl = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau#

import ipywidgets as widgets
from IPython.display import display, clear_output
from ipywidgets import HBox, VBox

# Widgets pour les constantes élastiques en GPa
mu_input = widgets.FloatText(value=80, description="μ = ", step=1)
mu_label = widgets.Label(value="GPa")

lambda_input = widgets.FloatText(value=120, description="λ = ", step=1)
lambda_label = widgets.Label(value="GPa")

# Organisation avec les labels d'unité
mu_box = HBox([mu_input, mu_label])
lambda_box = HBox([lambda_input, lambda_label])

# Fonction pour obtenir les constantes élastiques en Pa
def display_elastic_constants(change):
    mu = mu_input.value * 1e9  # Conversion en Pa
    lambda_ = lambda_input.value * 1e9  # Conversion en Pa
    
    # Afficher les valeurs en Pa
    clear_output(wait=True)  # Clear previous output
    
    return mu, lambda_  # Retourner les valeurs en Pa

# Liaison de la fonction d'affichage aux événements de changement de valeur des widgets
mu_input.observe(display_elastic_constants, names='value')
lambda_input.observe(display_elastic_constants, names='value')

# Affichage des widgets
ui_elastic = VBox([mu_box, lambda_box])
display(ui_elastic)

mu, lambda_ = display_elastic_constants('value')
print(f"Module de cisaillement (μ) en Pa: {mu:.2f}")
print(f"Premier paramètre de Lamé (λ) en Pa: {lambda_:.2f}")

Module de cisaillement (μ) en Pa: 80000000000.00
Premier paramètre de Lamé (λ) en Pa: 120000000000.00

Définition des chargements#

import ipywidgets as widgets
from IPython.display import display
import numpy as np

def create_angle_widget():
    # Création du widget curseur
    angle_slider = widgets.FloatSlider(
        value=20.0,  # Valeur initiale à 20 degrés
        min=0,
        max=360,
        step=1,
        description='Angle (°):',
        disabled=False,
        continuous_update=True,
        orientation='horizontal',
        readout=True,
        readout_format='.0f',
    )

    # Widget de sortie pour afficher les valeurs
    output = widgets.Output()

    # Fonction pour mettre à jour l'affichage
    def update_angle(change):
        ang_degrees = change['new']
        ang_radians = np.radians(ang_degrees)
        
        with output:
            output.clear_output(wait=True)
            print(f"Angle : {ang_degrees:.0f}° ({ang_radians:.2f} radians)")
        
        return ang_radians

    # Attacher la fonction de mise à jour au widget
    angle_slider.observe(update_angle, names='value')

    # Afficher le widget et la sortie
    display(widgets.VBox([angle_slider, output]))

    # Initialiser l'affichage
    update_angle({'new': angle_slider.value})

    return angle_slider

# Créer et afficher le widget
angle_widget = create_angle_widget()

# Fonction pour obtenir la valeur actuelle en radians
def get_angle_radians():
    return np.radians(angle_widget.value)

Taux de torsion (α) = 0.3491 rad/m
Rayon équivalent (R) = 0.2000 m
Rigidité en torsion (C) = 7.0184e+07 N·m²

# Calcul de la constante A pour le champ de torsion
A = 2 * C / (pi * R**4)  # A est un scalaire constant positif
print(f"Densité de contrainte de torsion (A_τ) = {A:.4e} N/m³")

# Définition des coordonnées spatiales
x = ufl.SpatialCoordinate(domain)

# Définition des composantes du champ de torsion (vecteur des Contraintes de Cauchy)
T1 = -A * x[1]  # Composante x du champ de torsion / Contrainte de cisaillement dans la direction x
T2 = A * x[0]   # Composante y du champ de torsion / Contrainte de cisaillement dans la direction y
T3 = ufl.as_ufl(0.0)  # Composante z nulle / Pas de contrainte dans la direction z (axe de torsion)

# Combinaison en un champ vectoriel de torsion
T = ufl.as_vector([T1, T2, T3])
print("Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)")

Densité de contrainte de torsion (A_τ) = 2.7925e+10 N/m³
Vecteur de contrainte de Cauchy défini pour la torsion (N/m²)

Résolution#

# Créez un widget Dropdown pour la sélection de la condition aux limites
bc_selector = widgets.Dropdown(
    options=[
        ('Aucune', 'none'),
        ('Face inférieure S0', 'S0'),
        ('Face supérieure Sh', 'Sh'),
        ('Surface latérale Sl', 'Sl')
    ],
    value='S0',  # Valeur par défaut
    description='Encastrement:',
    disabled=False,
)

# Dictionnaire pour mapper les sélections aux conditions aux limites
bc_dict = {
    'none': [],
    'S0': [bc0_S0],
    'Sh': [bc0_Sh],
    'Sl': [bc0_Sl]
}

# Créez un widget de sortie pour afficher les messages
output = widgets.Output()

# Fonction pour mettre à jour le problème
def update_problem(change):
    global problem
    selected_bc = change['new']
    problem = LinearProblem(a, L, bcs=bc_dict[selected_bc], 
                            petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
    
    # Mise à jour de l'affichage
    with output:
        output.clear_output(wait=True)
        print(f"Encastrement appliquée sur {selected_bc}")
            
# Attacher la fonction de mise à jour au widget
bc_selector.observe(update_problem, names='value')

# Afficher le widget
display(bc_selector, output)

# Initialisation du problème avec la valeur par défaut
update_problem({'new': bc_selector.value})

uh = problem.solve()

Visualisation des résultats#

Visualisation des déplacements#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats#

Calcul des valeurs maximales#

Moment d'inertie polaire J : 2.513274e-03 m^4
Angle de torsion théorique θ : 0.349066 radians
Angle de torsion théorique θ : 20.00 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 0.347123 radians
Rotation Z maximale simulée : 19.89 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 0.347342 radians
Norme de rotation maximale : 19.90 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 2.18 degrés
RX max  : 0.037984 radians
RY max : 2.13 degrés
RY max  : 0.037210 radians

Pb 2

Contents

Pb 2#

Problème torsion d’arbres#

Librairies#

Géométrie#

Définition de la géométrie et du maillage#

Configuration du problème#

Définition de l’espace fonctionnel#

Définition des frontière du domaine#

Visualisation des frontière du domaine#

Conditions aux limites#

Propriétés matériau#

Définition des chargements#

Résolution#

Visualisation des résultats#

Visualisation des déplacements#

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

Visualisation des contraintes de von Mises#

Analyses des résultats#

Calcul des valeurs maximales#